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ABSTRACT 

An intensive study of the assumed variable 
distribution necessary for the Assumed Displacement 
Formulation, the Hellinger-Reissner Formulation, and the 
Hu-Washizu Formulation is made in a unified manner. With 
emphasis on physical explanation, a systematic method for 
the Hybrid Stress element construction is outlined. The 
numerical examples employ four and eight node plane stress 
elements and eight and twenty node solid elements. 
Computation cost study indicates that the hybrid stress 
element derived using recently developed Uncoupled Stress 
Formulation is comparable in CPU time to the Assumed 
Displacement element. Overall, main emphasis is placed on 
providing a broader understanding of the Hybrid Stress 
Formul a tion . 
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1 . INTRODUCTION 


The Finite Element Method is a numerical method with 
firmly established mathematical foundation. Popularized 
by the broad applicability, large finite element codes 
play a dominant role in the current structural analysis. 
Thus, in the past years, intense effort was applied to 
improve and optimize the finite element method* 

The main thrust of this report is to present an 
explanation of the finite element method in the view of 
the theory coupled with computational algorythm* The 
functionals considered are that of the Hu-Washizu 
principle , TThv' the Hellinger-Reissner pr incipl e , 7Tf\ , and 
the principle of minimum potential energy, 7Tp • By a 
direct comparison of these three functional, the role of 
the assumed field variables can be clarified* 
Furthermore, through a comparative evaluation, these 
theories will be approached in a unified manner* 

The motivation came from earlier attempts for the 
explanation of the class of elements derived under hybrid 
stress method. Several months of numerical and literature 
research demanded more systematic method for the research 
into the hybrid stress elements* Using the latest 
development in the hybrid element research, a systematic 
method will be constructed step by step starting from the 
governing equations of elasticity* 
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2 . GOVERNING EQUATIONS OF ELASTICITY 


Since finite element methods in solid continuum 
operate on the governing equations of elasticity, several 
key observations must be emphasized before introducing the 
weak or the variational form of the equations. In this 
analysis, only the small displacement theory of elasticity 
will be considered. Also, the Rectangular Cartesian 
Coordinates will be employed for defining the three 
dimensional space. 


STRESS 


The condition for the stresses are obtained directly 
from the Newton's laws of motion pertaining to the bodies 
at rest, i.e. the force and moment equilibrium. 

Application of these laws on the stress yields the 
Equations of Equilibrium for stress. 
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where , 

U = Six stress components * Gy ®xy ^y£ 
F * Body force components 
D * Matrix of differential operators 
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Thus 

, the chosen 

stress 

field must 

satisfy 

the 

pointwise 

homogeneous e 

quil ibrium 

dictated by 

the Newton' s 

laws* This applies to 

the any 

functional 

whether 

it 

contains 

stress as 

expl ic it 

unknown 

functions 

or 


implicitly through the unknown displacement functions* 
Further discussion is reserved for later sections* 


DISPLACEMENT 

Since the class of elements presented require 
continuity, the isoparametric formulation is the obvious 
choice* These elements will be later refered to as the C° 
continuity elements* The procedure in obtaining the 
interpolation functions are well outlined in most finite 
element reference books, e.g. [2]* The notation used for 
the displacements are 

( 2 . 2 ) 



STRAINS 

The strains are defined as 
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These six independent partial differential equations will 

be denoted as the Strain-Displacement Relations. 

For the functionals, 7T P and ^ Pk. • where 

displacements are included explicitly as unknowns in 

formulating the elasticity boundary value problem, no 

further conditions are necessary to ensure the existence 

of single-valued displacement. However, for 7T U¥Ai / since 

n w 

the strains are also considered as independent unknowns, 
the compatibility conditions must be imposed on the 
strain-displacement relation to ensure that the unknown 
strains are indeed compatible with the unknown 
displacements • 
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These equations are known as the St • -Venant 1 s 
compatibility equations. However the six equations do not 
represent six independent conditions. Yet the usual 
procedure is to include all six equations in the interior 
problem formulation, but to remember that they represent 
only three independent conditions [3]. 

Thus far, all the equations presented are completely 
independent of the relationship between stress and strain. 
They are applicable to any type of continuous body 
undergoing small displacement. However, to predict the 
behavior of a structure it is also necessary to know the 
components of stress as functions of the components of 
strain and vice versa. Through the Stress-Strain Relation 
the material properties of the body enter the problem. In 
the following development, the materials considered will 
be assumed to follow the generalized Hooke's law with 21 
independent constants. 

The stress-strain relation can be written in a matrix 
form as 

<T = C £ 


and 


! = S Z 


where , 



( 2 . 5 ) 


For this stage the most predominating factor must be 
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pointed out. This is simply that the stresses, the 
displacements# and the strains must satisfy these 
governing equations. When the approximate functions are 
used to solve these equations, the optimum choice of the 
assumed function is the one that satisfy each of these 
equations a priori to the highest polynomial order. This 
argument can be used to sort through all the admissible 
functions in formulating the best general finite element. 
Thus# whatever variational functional used to formulate 
the element# one must always refer back to the physics of 
the problem and not sink into the mathematical generality. 


n 
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3. VARIATIONAL FUNCTIONALS 


Without belaboring on the actual derivation of each 
functionals, the main references will be indicated and 
only the pertinent information for this discussion will be 
outlined. Presentation is organized to clearly indicate 
the a priori conditions and the stationary conditions 
which provide the Euler equation for each functional* 


PRINCIPLE OF MINIMUM POTENTIAL ENERGY , TTp ( u) 


A priori conditions: 


£ *.D u. 
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ic - a 


} 


in 


V 


on 


S u 


(3.1) 



e T c t n 

«v - 


( £ T S «•* - ( J r £ 


(3.2) 


The 


stationary condition, 
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where , 

€ = Strain 

0" = Stress 

C = Stress-strain law 
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U = Displacement 
IX = Prescribed displacements 
T = Prescribed traction 
F = Body force 

V =* Volume 

5^ - Portion of surface where displacements 
are prescribed 

C = Portion of surface where stresses 
<? 

are prescribed 


In 7T P , the strain-displacement relation, 
stress-strain relation, and prescribed displacements are 
identically satisfied. However, the equilibrium condition 
and the prescribed tractions are only satisfied in a 
variational sense. By relaxing all three a priori 
conditions by using the Lagrange multiplier method, the 
generalized functional, 7T W ^ , can be obtained. 

HP— WASHIZU PRINCIPLE , <£'£'£•) 


A priori conditions: NONE 


7T h „ = ( [i e r i£ - r«l dv - \ x T u ds 
- ( s’Cs-es) ^ - ( z 

J y J s u 

The stationary condition, <£tT hm/ =0 , yields 
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U =. LL on 5 U 

-V *v> 

I = 'itz = T on 5, 

®* ■* ~^y '*’ 

•i K = ••• D>rec4i»nai Cosines 


All governing equations are satisfied in a 
variational sense. Yet the price paid for the 

generalization is the two additional unknowns € and 
<T . The Hellinger-Reissner Principle can be easily 
obtained by introducing the stress-strain relation* 


HELLINGER-REISSNER PRINCIPLE , 7Tj^ (u, £ ) 
A. priori conditions: 



(3. 6) 
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The stationary condition, <£7T^ =0 , yields 
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T = T on 

U =. U on 5 U 


By satisfying the equilibrium equations identically, 
7T r and can be shown to reduce to the Hybrid 

formulation via Principle of Complementary Energy [4] • 
Yet with availability of isoparametric formulation for C Q 
continuity elements, the Hybrid elements can be formulated 
using TT^ or ^hw in wore expedient manner* 

In a recent publication by Pian and Chen [5] , new 
formulation for hybrid elements provides a method to 
directly take advantage of the sparse nature of matrices 
involved in the construction of the hybrid elements* 
Briefly, the new formulation introduces an additional 
displacement field which act as a Lagrange multiplier on 
the homogeneous equilibrium equations* The major 
advantages will be described in the later sections. These 
new functionals will be refered to as the uncoupled stress 
versions • 

UNCOUPLED STRESS VERSION OF 7^ AND 7T HW 


The element displacement u is divided into two 
separate parts. 


“• * Ht + Xx 


(3.9) 


u g is the usual compatible displacement in terms of the 
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nodal displacements • u ^ is the additional internal 
displacement which acts as the Lagrange multiplier for the 
homogeneous equilibrium* When this form of displacement 
is implemented into the and TT Hw , Equations (3.7) and 
( 3 . 4 ) , with body force and prescribed traction terms 
dropped, the resulting equations are 


tt r = ( [-i2 T s»- * -(s. T 2) T “J d ' / ( 3 , ®> 

fji £ t ££ - 2 t £ * 2 t C5“j) * S-Vdu,)] dv (3.ii) 


As a remark on notation, since the hybrid formulation can 
be obtained via or TT HW *- n both coupled and uncoupled 
stress versions , ^ will be used to refer to the general 
hybrid element case where the stresses are constrained to 
be in equilibrium. 
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4 . FIELD VARIABLES and STATEMENT of EQUIVALENCE 


The first question that must be answered in any 
finite element development is whether the method 

converges* The mathematical proof for the convergence of 
7T P and are given in the References [6] and [7]. 

Yet these proofs only justify the use of the finite 

element method and not whether it is feasible for 

engineering application* The demanding requirement is the 
rate of convergence* 

The two modes of convergence are the h-convergence , 
the diameter of the largest element, h-max , approaches 
zero, and the p-converg ence , the minumum order of the 
polynomial basis functions, p-min, approaches infinity. 
The comparison of the rate of convergence between these 
two modes for -TTp is covered by Bubuska and Szabo [8] • 
The results given in this article should be clarified in 
the view of general purpose elements* First, most 

structural problems only require low order polynomial* 
Thus the lower order elements can be employed for most 
practical engineering problems* Secondly the special 
purpose elements and other methods such as the reduced 
integration techniques should be considered. The main 

reason for the discussion of this article is to point out 
that before the general purpose elements can be 

constructed, better unified understanding of the special 
purpose elements and various other techniques is 
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necessary* As the solution of the problem by finite 
element method becomes more complex, as exemplified in 
crack and composite material analysis, the analyst 
judgement and experience will not be adequate in 
determining the accuracy of the solution. 

Another important point that must be defined involves 
the past attempts to improve elements based upon / Tp 
The first workable attempt was presented by Wilson [9] 
with his incompatible displacement models. Recognizing 
the deficiency in the isoparametric displacement elements, 
he introduced the bending modes into the element 
displacement field. Other trials that falls into this 
category can be exemplified by the introduction of 
singularity by distorting the nodes. The second category 
can be classified into the scheme of reduced integration 
technique. This technique works well to prevent locking 
in the generalized shell element as well as other 
applications. The major point that most researchers in 
finite element method ignored is the fact that both 
Wilson's incompatible element and the reduced integration 
technique result in identical stiffness matrix as the 
hybrid formulation. This bold clue that hybrid 

formulation can reveal the way to develop the optimal 
element has been casually dismissed in view of 
computational cost in a premature fashion. 

With this emphasis clearly made, a close examination 
of the finite element method in a unified manner can be 
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made. First the field variable for each functional will 
be listed. 

t displacement 

7T^ displacement, stress 

7T^ (“,£,0 displacement, stress, strain 

The Statement of Equivalence describes that when the 

displacement modes from TTp yields the stress and strain 
modes for TT^ and TT HW , the three functionals are 
equivalent or identical. To illustrate, take the simplest 
case of linear displacement modes. Since this yields 
constant stress and strain modes, if <T and £ contains 
the constant modes, the three functionals are equivalent 
up to a constant stress and strain. The argument directly 
leads to the order of equivalence. 

From the consequence of above statement, using the 
inductive reasoning, any class of problems that have 
constant stress of strain will yield the exact solution 
using any of the above functionals. Using the Statement 
of Equivalence as a benchmark, the examination of each 
functionals separately may be continued. As a warning, 
the mechanics of construction of the desirable field 
variables are left for the later section and should not 
enter here as a consideration. This is justifiable since 
for clarification purposes, ideal cases may be presented. 
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~ TT P (u) — 


A priori conditions (3.1) imply that the choice of u 
automatically determines Z and <£ . Under isoparametric 
formulation , there is no flexibility on the choice of u. 
Therefore, further discussion on u is not necessary. 

Recall, from the section 2, that the governing 
equations are the equilibrium equations (2.1), the 
strain- displacement relation (2.3), and the stress-strain 
relation (2.5). Eq ua tions (2.3) and (2.5) are 
automatically satisfied under the isoparametric 
formulation. However, the equilibrium equations (2.1) 
applied to the stress components obtained from the 
displacement modes is not satisfied. Furthermore, the 
stress modes are artificially coupled. Artificial in a 
sense that one stress component coupling into another in a 
way not possible under equilibrium considerations. The 
undesirable trait appears as locking for the bending 
problem. The solution to bypass this problem is to assume 
separate u and £ with Z satisfying equilibrium, thereby 
satisfying all of the governing equations. 

— TT_ ( u , 2 ) — 


Under the governing equations the forces , • an ^ the 
deflections, u, couple only through the stress-strain 
relations. By choosing u and independently, all three 
governing equations can be satisfied a priori. Therefore, 
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with u obtained through isoparametric formulation and J? 
chosen correctly, an optimal general finite element can be 
con str uc ted • 

In selecting the stress field, two articles provide 
means to overcome the preliminary roadblocks. A method to 
bypass the zero-energy deformation modes due to the rank 
deficiency in the stiffness matrix is outlined by Plan and 
ChenflO]. Their method simply matches a stress mode for 
each possible strain mode in order to force non-zero 
strain energy. This method is easy to apply and provide 
effective means to detect and eliminate the troublesome 
zero-energy deformation modes. 

The second article, by Tong and Plan [7], on the 
convergence of the finite element method based on assumed 
stress and displacement distribution, outlines the 
procedure to gage the accuracy of such an element. In 
order to obtain progressively better accuracy, both the 
stress and the displacement approximations must be 
improved properly and simultaneously. In another words, 
the largest error, whether the error is from the u or (J* , 
is the error of the finite element approximation. 

The Hu-Washizu principle is a generalized functional 
which allows full flexibility on the choice of the assumed 
field variables. This is called "generalized" since no a 
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priori assumption is made and thus the governing equations 
of elasticity are satisfied as Euler equations of the 
functional. Even with the full flexibility# the choice of 
the assumed field variables follow the same constraints as 
indicated for 7T^ . Although using to obtain a 
hybrid element seem a "round-about" way, the reason will 
be apparent under computational efficiency dealing with 
anisotropic material. 

ROLE OF THE LAGRANGE MULTIPLIER 

In the finite element method, the Lagrange multiplier 
relaxes the corresponding governing equation. Thus, the 
resulting functional becomes more flexible for the 
implementation of the desired element properties. From 
the view of understanding the mechanics of the finite 
element method, the Lagrange multiplier decouples the role 
of the assumed variables. From the view of computational 
algorythm, the Lagrange multiplier introduces flexibility 
needed to reduce computational cost. By decoupling the 
assumed variables, algorythms can be implemented to take 
full advantage of the sparse nature of the necessary 
matrices for element generation. 

Even with added flexibility, emphasis must be made on 
a simple knowledge that for a given well posed boundary 
condition, the governing equations of elasticity has a 
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unique solution. Thus restriction depends on whether the 
element is general purpose or special purpose element. 
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5. GENERAL PURPOSE FINITE ELEMENTS 


The key to constructing the general purpose finite 
element is the actual understanding of each element in 
both mathematical and physical sense* In the early stages 
of hybrid element research. Irons provided a strong 
physical insight in the process of the development of an 
assumed stress version of the Wilson's incompatible 8-node 
isoparametric brick element [11] • He obtained the correct 
stress parameters by simply choosing the modes that 
describes physical states of the classic problems* For 

example, the pure bending modes should be included in the 
stress field, whereas other terms that contribute to the 
spurious strain energy should be excluded* As Irons 
pointed out, the isoparametric element prevents the need 
for engineering insight* Both the researchers and 

analysts must recognize the limitations of each element 
which is clearly provided in the development process of 
hybrid/mixed formulations* The closed minded approach of 
using the simplest functional, iTp , and ignoring the rest 
will hinder the progress. Recognize that the construction 
of isoparametric element is simple because the steps allow 
no flexibility and thus faced with its full limitations. 

Falling back on the Statement of Equivalence, any 
desired element characteristics, when proven to exist in 
TTp formulation can be reproduced in the other 
functionals. Since the Wilson's incompatible element and 
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the reduced integration scheme can produce pure bending 
modes, these desired characteristics can be reproduced 
using T T R or TT^ W . The difference lies in the algorythm 
used for the construction of each element. 

Thus far since only the C° continuity elements were 
considered, the time is ripe for the discussion of the 
elements requiring C l continuity as later referred to as 
C 1 continuity elements. Since the original application of 
the hybrid formulation was intended for the C* continuity 
elements using the modified complementary energy principle 
[12], this discussion inevitably follows. 

First recall that the primary reason for the 
derivation of the beam, plate, and shell theories were to 
simplify the analysis in order to obtain an analytical 
solution. Each theory is based upon the assumption that 
one or more dimensions of the problem collapses. To 
illustrate, the governing equations can be 
non-dimensionali zed using a characteristic dimensions of 
the problem. As one of the dimensions collapses, for 
example the thickness, the equations can be expanded into 
a perturbation series. Thus, the structural theories are 
basically the leading order, or the zeroeth order, 
equation of the governing partial differential equations. 
As the perturbation parameter increases, the accuracy of 
the leading order equation diminishes. In order to 
improve the range of validity of the approximation, later 
works introduced the transverse shear effect as 
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exemplified by the Mindlin plate theory. Although the 
attempt will not be made here, the transverse shear effect 
probably is the first order correction of the perturbation 
expansion • 

The price paid for the simplification made by the 
structural theory is the requirement for C* continuity in 
the finite element analysis. Even to this day the 

agreement whether this price is justified in the finite 
element analysis has not been reached. However, the 
degenerated plate and shell element is gaining popularity 
for the practical applications. Note that the degenerated 
element class only require C° continuity for the 

displacements . 

In the article by Bathe and Bolourchi [13] , an 
extensive coverage of the degenerated plate and shell 
element using "JTp formulation has been made. The 

numerical examples given indicated that the best solution 
is obtained by using the selective, or reduced, 
integration technique. Thus by using the hybrid 

formulation the results obtained through reduced 
integration can be identically reproduced. 

The major advantage of the C 1 continuity elements is 
the ability to represent the bending behavior. If the C° 
continuity element can represent the bending behavior, 
then the general purpose finite element can indeed be 
constructed. The line of research along this path using 
7T p formulation has been hindered through the difficulty 


21 



arising with the locking problem. Since locking is caused 
by artificial coupling that is inherent in the 
isoparametric assumed displacement formulation, the 
problem can be easily remedied by the use of the hybrid 
formulation. Furthermore, the bending modes can be 
implemented even in the linear hybrid elements. Then the 
only limitation on using a solid element to model, for 
example, the thin plate behavior is the numerical 
stability. By using high enough decimal precision this 
limitation can be avoided. 
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6. SYMMETRY AND RELAXATION 


As numerically demonstrated by Plan, Chen, and Kang 
[14] , the symmetry condition is an important criterion 
that must be considered in obtaining the assumed field 
variable distribution. This is physically consistant 
since the finite element should be symmetric in all three 
directions. For example, when choosing the stress modes, 
bending behavior in all three coordinates should be 
represented. This criterion is invaluable tool for the 
solid element construction. 

Before the relaxation condition can be described, a 
better understanding of the zero— energy deformation mode 
(ZEDM) is in order. Although mathematically ZEDM is a 
rank deficiency beyond the rigid body modes in the 
stiffness matrix and its prevention outlined in the 
reference [10], more physical explanation is necessary. 
In the case of solid elements, each node contains three 
independent displacements, u, v, and w, to describe all 
possible motions of the node. With all the element nodes 
moving in conjunction, all possible deflection modes can 
be determined. For most of these deflection modes 
physical significance can be attached such as pure 
tension, shear, or bending. For the rest, no such 
physical association can be made. Granted each material 
point has three degrees of freedom, but that point cannot 
be considered as an isolated particle in free space. A 
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continuum collection of material particles has additional 
physical deformation restrictions. The modes with no 
physical association arise strictly through the 
mathematical modeling process. 

Whether the deflection mode has a physical 
significance or not, the corresponding non-orthogonal 
stress mode must be provided in order to prevent ZEDM. 
Non-orthogonality of the modes guarantee non-zero strain 
energy. Therefore by identification of all the deflection 
modes, ZEDM 9 s can be easily eliminated. 

As a classic example of the duality principle, a 
parallel analogy is presented by Loikkanen [15] using the 
stress modes. In the article, he nicknamed the "nonsense" 
stresses refering to the stress modes with no physical 
association. However, for the purposes of generalization 
for the higher order elements, the use of deflection modes 
will be more convenient since through the isoparametric 
formulation the deflection modes are given and the stress 
modes yet to be determined. Note that any non-orthogonal 
stress mode can be used even though the lowest order mode 
is preferred for the numerical integration considerations. 
Figure 6.1 pictorially summarizes the above discussion. 
For continuum problems, a complex combination of forces 
required to excite such an element deflection within a 
mesh arises only for the occasion when the mesh is much 
too coa r se . 

From here on, the "nonsense” stress mode refers to 
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any stress mode with only general non— orthogonal ity 
restrictions used in conjunction with the deflection modes 
with no physical association# With this in mind the 
relaxation condition can be simply stated# For the 
"nonsense" stress mode and the corresponding deflection 
mode # the governing equations can be relaxed without loss 
of accuracy of the element# In the practical application 
of the above statement# in M ^ the equilibrium condition 
can be relaxed for the "nonsense" stress modes and in TT HW 
the stress-strain relation can also be relaxed# 
Surprisingly large computational costs can be saved by 
using the relaxation condition# In addition# more 
accurate solution can be achieved by eliminating all the 
supporting terms necessary to keep these "nonsense" stress 
modes in equilibrium* 
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7 . FINITE ELEMENT METHOD 

For the finite element method presented below, only a 
single element domain is under consideration for each of 
the functionals. 

7.1 PRINCIPLE OF MINIMUM POTENTIAL ENERGY , 7T p ( u) 

In using the principle of minimum potential energy to 
formulate the element stiffness matrix the following 
functional for an element should be stationary, 

TTp - i f - \ Jt T H dv - [ X T ~ dS C3.2) 

r <v J v 

In the matrix form, 

-ir. * i * T *t - « T * tr-o 

where 

K =* stiffness matrix 
% =■ Nodal displacement 

Q - Load matrix 

Furthermore, 

S 5 5, £ T £ § ^ (7.i) 

q - S a r Sdv ♦ ( * T ! ^ 

^ a 

The relations used above are 
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U = A/ % 


( 7 . 3 ) 


£ = DU r (D A /) Z = B f 

-s ' — ** ^ ~ «v -w 

where 

Ai = Interpolation matrix 

D = Derivative matrix 

fl = Strain-displacement matrix 
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7 . 2 HELLINGER-REISSNER PRINCIPLE , 7 T_ ( u , £ ) 


For the Hell inger-Reissner principle , the following 
functional for an element should be stationary, 


*tt = ( [--k + e T c^if) -£ T «] m 

* 

. j fujs - ^ X T { ) As 


C 3 . 7 ) 


In the matrix form, assuming u^u on S 4 


JL 


where 


and 


/£ T "£ -/£ t S£ ' 


H = 
Q = 
<? 


cr = p 
u = N % 

( P T 5 P d V 
f P t 5 JIT 

' iy 


= f /V T F civ + f ^ T I 
j v ~ J S„ 


T As 


C 7 . 4 .) 


(7.0 


(7.0 


The matrix, P, contains the internal stress modes and 
jS is the corresponding unknown coefficients. In order to 
obtain the standard stiffness matrix form, use the 
stationary condition on Equation (7.4) with respect to ^ • 
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-*-rr K 


_f> H 


fi - 5 f 


= o 


(7.7) 


or 


/ = H' <7 ? 

yw ^ ^ 


Substituting back into Equation (7.4) yields the familiar 
form 

I - T ~ — T 


* -£ t K 9 - Q ? 

^ -s. -W 


(7-?) 


where , 

K = o T h h q 


(7.9) 


In order to obtain the stiffness matrix, Equation 
(7.9), H matrix must be inverted. The order of II matrix 
is the number of stress parameters, ^ * s, used for the 
element development. The condition for the existence of 
the solution for ■ s is given in the reference [7]. Let 
m equal to the number of ^'s, k equal to the number of 
element degrees of freedom, and 1 equal to the number of 
rigid body modes. Then this condition is simply that 
m fct k-1 • Therefore, since the order of the H, matrix is 
(m«m), the minimum possible m is k-1. Since the number of 
algebraic steps required for the inversion is on the order 
of m* , the computational cost for generating the stiffness 
matrix can be significantly more than that for the 
element stiffness matrix. 
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7 . 3 HU-WASHIZU PRINCIPLE , TT HW ( U , £ , £ ) 


For the Hu-Washizu principle, following 
~fT HW should be stationary. 


TT 


HW 


= J [x e T c£ - F ] J / - J T T u ds 


- { £ T (£ ' 5“ ) ^ " f I T ^ - 5 

S u 

In the matrix form, assuming u=u on , 

t r = 4 - 8 r H oc + f 

* ' Hw *’ '" w /^v ^ *w ^ *v 


*) «= 


where 


s* = r* 

£ = p* 

u = £/ f 


and 


H = 


7 

-V 


<5 

Q 


f P T P dv' 

Jy ^ ^ 

=. J p r c p jy 

= ( p t b jv 

Jy 

= j /V T FdV + j" /f ds 


f unctional 


(3. + ) 


(7. /o) 


C7.») 


(7.'Z) 
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With the assumpt 
the generalized 
both stress and 
linear combinat 
major deficiency 
that both the 
stress- strain re 
physically signi 
Using the s 

& > 

v 

or 

<* - 


£ 


Substituting bac 


^HW 


X 

z 


where 


ion that the material behaves according to 
Hooke's law, the use of same modes for 
strain is consistent since strain is a 
ion of stress and vice versa. However the 
in using TT HW formulation in this form is 
stress equilibrium condition and the 
lation cannot be enforced a priori for the 
f icant modes • 

tationary condition with respect to oc and 


T ck - hlfi - O 
-Hoc -+ c,t = O 

- -w 


C 7. > 3 ) 


H" 1 Z - H 1 T H ' <5 % 


(7. 14) 


k into Equation (7.10) yields 

l * i “ § T I C 7 '5) 


K = H 1 J ri Cf 


(7. U ) 
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7.4 UNCOUPLED STRESS VERSION OF TT R AND TT rtw 

As given by Plan and Chen(5] uncoupled stress 
Of 7TV with additional Lagrange multipliers u ^ is 

t r R = \ [-i<l T s£ + ~ ( a* 1 ^ 

In the matrix form, 

t* = ' i /? T «/£ * £ t 5I ' /£ T « i 

where 

£ - £ £ 
if = + «a 

ff* = 

“* * " A 

D t c; ~ E jff 

D T = Homogeneous equilibrium operator 

and 

H =■ \ P T -S P dv^ 

Gf - \ P T B d V 

A = 5 y £ t * dv 


(3.10) 


( 7. • 7 ) 


[7.1S) 


(7.H) 
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Using the stationary condition with respect to ^ and , 


= - R £ ■+ £} f - R "X =o 

-*{3 ~ ' ~ ~ " " 

r - R T £ = O 

3 A - ^ 

or 

/€ «" (§1 - *£) 

* = (b t h"p.)‘ fCn^Cr t 


Combining , 


a 

where 

Q 

mm 



$ t 

<*%* -v 


§ - 6(e T ti"d)'’s r ^"9 


The resulting stiffness matrix is 

K = § T H ^ 


C7. 24>) 
(.7- 2 1) 

C 7 - 22 ) 
(7. 2 B) 

(7. 2 4) 


Although at a first glance the uncoupled stress 
version of / seems overly complex, this functional 
elegantly takes a full advantage of the sparse nature of 
each matrices. Observe that the necessary bulk matrix 
multiplication. Equation (7.23), is performed only once. 
With an appropriate choice of , this formulation 
becomes identical to the standard 7T^ , Equation (3.7). 
More careful examination of the mechanics of the element 
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construction is worthwhile at this point* 

Since the equilibrium condition is imposed through 
the use of the Lagrange multiplier u^ , the stress 
components are uncoupled. 


V 


p, 


a 

h 



^ o 

Iz 


b 


r — 

« ** 

i 

t* 

is 

* - 


Zc _ 


w - 


(7.25) 


For the reasons of computational efficiency, / 1ft should 

be only used for the isotropic material. For the 
anisotropic material , the uncoupled 7Tj!| W will be much more 
efficient . 

From the Equations (7.24) and (7.23), to generate the 
stiffness matrix H and (R^H^R) matrices must be inverted. 
First examine the H matrix for the isotropic material. 
Define 

P; = ( (7.24) 

J v 

After algebraic manipulation, from Equation (7.19) 
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H-= E 


1 -0*- 


l"' 


assuming 


- -I 

-i) P. 


f, 


Sym. 


* f. 

f," 


sii p" 

2 


- p -. 

Z - 


L=* p“ 
2 


r. s r* s Pi 


( 7 . »7) 


Thus the order of inversion is drastically reduced* The 
choice of setting £,=£*=£3 instead of equating all six 
is induced by the necessary inversion of ( R r H R) matrix. 
The order of { R 1 " R) matrix is determined by the required 
number of A's, Equation (7.18), to impose equilibrium 
condition. The number of ^*s required can be determined 
by imposing equilibrium on the stress components and 
counting the number of constraints on the ^*s. Therefore 
the tradeoff is that if all six 1 s are equated, the 

number of ^*s will increase drastically. 

The reduction of the inversion order of JH matrix 
solves the question of the high computational cost using 
hybrid elements. An additional flexibility in the 
uncoupled stress formulation provides more versatile 
element. When the element is skewed the displacement 
modes becomes correspondingly skewed in the rectangular 
Cartesian coordinates. In order to exactly match these 
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skewed mode s , 


the assumed stress modes should also be 


defined in the natural coodinate system* The uncoupled 
stress formulation provide this flexibility. 

In order to establish the mechanics to assume the 
stress modes in the natural coordinates, further 
explanation of the R matrix, Equation (7.19), is needed. 
The stationary condition. Equation (7.20), and the 
Equation (7,19) provide two possible ways to obtain the R 
matrix • 

5 - L £ t " ** C7..9) 

= ° (7. ZO) 


Using Equation (7.20) the constraints on the f %' s obtained 
through the homogeneous equilibrium equations can be 
directly imposed if the stresses are defined in the 
rectangular Cartesian coordinates. The R matrix in this 
case simply zero or couple the appropriate f3 ' s. By 
assuming £ and u^ modes in natural coordinates and using 
Equation (7.19), similar type of constraints can be 
imposed • 


For the Hu-Washizu formulation. 


~JT hw = (" [x e T cz - £ t £ - o- T (^y) + M O n 

1 Jv 
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In the matrix form, 


ir H „ - * i'zs -fCBi '£st - 

where 

Z * if. 

£ o Pol 

u * «j <• “» 

?t‘ £>* 

ii » = M A 

^ /» -N- *s- 

5 T f = ?/- 

and 

P? = f P T P of^ 

v ~ ~ 

i = r p t c p d/ 

__ \ *V. -**■ <»■< 

^ J V 

(7 = f P T S dl/ 

— *v ^ 

= ( E t /M cj V 


C 7 . ) 


£ 7 . 2 ?) 


( 7 . 3 0 ) 


Using the stationary condition, 

- HC<- + (7 f ~ R% ~ 0 

T~ -- ^ •» *>. -v- 

lZEa. u = To< - h & - a 

~~ «w f^ 

3 a. 


(7.31) 
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o 


or 


3 2 s 



* - H ’’ [ § I - « ^ ] 

J9 = C 7 5 i> 

^ = (R T H H T B“g) 'g T H ■* J 


Define , 

W = H '* 3 H "’ (7. 3J) 

g » g - « ( 5. T iy2 )' " i? T jy £ (734} 


The resulting stiffness matrix is 

K = § 


( 7 . 3 S-) 


Since the Hu-Washizu formulation relaxes the 
stress-strain relation, further constraints on the choice 
of the P matrix must be established. These constraints 
arise due to the assumption that 


- “ f /I (7. 2?) 


Define P as 
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p ^ 


r. 


& 


5 




!■> 




( 7 . 3 *) 


For the anisotropic material law, in order to obtain 


<y = c e 


C z. s) 


exactly, the P matrix has to be chosen such that 


r» = tz ~ ~3 ” = ^ _ 

For the orthotropic case, only P 4 » P ^*2. 4 must be equal to 
each other. To satisfy the pointwise equilibrium, the 
procedure used for /f^ also applies to this functional# 

When the matrices required to generate the stiffness 
matrix for the uncoupled stress 7Tj^ and ~TT H „ , the 
equivalence condition for the two functionals are observed 

H (7.37) 


to be 


VJ - 


H ' 3 


_ -i 

H 


At this point, the purpose in exploring the 
principle can be easily demonstrated, 
significant difference between the fini 
formulation using TT^ and Tf^ arises for the 
material. To illustrate, for a 20-node solid 
minimum number of ' s is 54. In order to 


Hu-Washizu 
The most 
te element 
an i sotropic 
e 1 ement , the 
satisfy the 
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cubic ) 


requirement that all are equal, 120 ^ % s (full 

are required. To apply equilibrium and reduce this to 54 
^S's, 66 P\'s are necessary. Thus the order of (R t H~*R ) 
required for the stiffness matrix is 66* This defeats the 
purpose of using the uncoupled stress approach. 

In the Hu-Washizu formulation, the stress-strain 
relation can be relaxed for the "nonsense" stress terms 
analogous to relaxed equilibrium condition for 7T^ . 


Thus , 

the stress 

-strain 

parameters P. can be 
- 1 

assumed to b 

o nl y 

identical 

up to 

a chosen order • 

Furthermore 


different P^'s can be used since H is in a diagonal form. 


H 


-i 



(.7. 3S) 


where 


P-. 


— LL 


[ P T P. Jl/ 

J v - 1 


(c noi Summed) 


( 7 . 39 ) 


By also 

relaxing equilibrium for 

the "nonsen 

te rms , 

the total number of 

s required can 

reduced • 

A numerical example is 

provided in 

section. 




se" 
be g 
the 


stress 
rea tl y 
later 
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8 . ELEMENTS 


8 . 1 4-NODE LINEAR PLANE ELEMENTS 



ISOPARAMETRIC ASSUMED DISPLACEMENT ELEMENT 


Interpolation 
LL r 
- 


functions for 4-node 


w t 


plane elements are 

(f?o 

i‘ ',*. 3,4 


By expanding the interpolation 
displacement modes can be obtained* 
each mode is denoted OC, • 


functions the 

The coefficients of 


u = + °<zj + *3 7 + * 4-3 7 

v = * 5 + <* g ^ f + *8 37 


Without loss of generality, since isoparametric mapping 
has one-to-one correspondence, consider a rectangular 
element with x and y correspond to "J and *1 , 

respectively. Using the strain-displacement relation, the 
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strain modes corresponding to Equation (8.2) are 


** °<*y 


*7 + 


r xy » (*3 + * 6 } + + 


Isolating each strain modes, 


C, ' *2 


fy - *7 


El.ON(jAT/OA/ 

= *5 * ** } distortion 


(*• 3 ) 


( 9 . 4 ) 

( 9 . 5 ) 



£* - o 

ty = «, * (9.i) 

*x y s «r y 


The two modes indicated in Equation (8*6) demonstrate the 
artificial coupling inherent in the isoparametric 
formulation. This coupling nature is the reason why 
elements handles the pure bending poorly* In pure bending 
the shear term approaches zero as the thickness 
diminishes* The two modes in Equation (8.6) with the same 
coefficient must increment both normal and shear strains 
s imul ta neo usly • 


WILSON'S INCOMPATIBLE ELEMENT 


Wilson's incompatible displacement fields for 4-node 


plane element are 
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u - «, + «»5 f » "+ 31 + tj'f) 

V= * 5 .« <5 +«,7 ♦ «iil + V'** 1 + 

The four additional terms are the incompatible 
displacement terms to represent the bending behavior. The 
analysis of this element follows the same procedure. 

i , * ** + ot 4y - z *« * 

. w (*.*) 

€ y C <X 7 + « ? x - 2 \ y 

y x ^ •- («<jt ^ + + °<sy "^^y* 2 ^3 x 

By redefining the coefficients in Equation (8.8) as 

A * r ' 2 * 3 (*.<?) 

A z = °<s ~ 

the strain modes can be isolated. Upon sub stitution # 

€ x = + A,y + z?\ 2 y - z* t x 

£^ = <x^ +■ A* * " z ^4 / (F. JO) 

r xy * + + -A. 1 * 4 y 

Although the bending behavior is retained in the Wilson's 
incompatible element by the added Ajand modes, the 

element fails the patch test unless the shape of the 
element is rectang ular [ 1 1 ] . 
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HYBRID ELEMENT 


A demonstration of constructing the uncoupled stress 
TT^ and elements is left for more general solid 

element discussion* For the plane elements, only the 
hybrid elements constructed using the standard 7["^ are 
ill us tra ted • 

To eliminate the zero-energy deformation modes, the 
strain modes obtained from the displacement modes. 
Equation (8.3), provides a convenient method. Each strain 
mode must be matched with stress mode to form non-zero 
energy. The technique is to first eliminate artificial 
coupling by using additional ^5 1 s in the matching process. 
Thus Equation (8.3) forms 

- /*. r A y 

°y * fix * Ps * 

-- h * P‘ x * fr y 

Other modes, such as OJ *x and Oy «y, are not 

considered since they introduce additional error. Since 
these modes are not present in the strain modes, they 
interact with all the non-orthogonal terms in the element 
energy integral and thus introducing artificial coupling 
error. To illustrate, if the mode 0^ * x inc l ude<i 

equilibrium condition requires that ^xy * -h y . The 

coupling arise when <J^y = -^y interact with mode in 
Equation (8.3). Thus the shear again couple into the 

bending behavior. 


(S.,1) 
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Applying the equilibrium condition on Equation (8*11) 


yields the optimal assumed stress modes* 

f - 7 * ° 

/*< * ° 

* p, * Y 

®v = A *■ * 

°"«y ’ /* 


(?•'*) 


Since these stress modes are not invariant with respect to 
the reference coordinates, local axes should be used* 
Forcing invariance through the use of complete 
polynomial [ 16] diminishes the true advantage of the hybrid 
formulation for the general purpose elements. 

At this stage the reason why the selective 
integration technique and the Wilson’s incompatible 
element reduce to hybrid element for the rectangular 
geometry is clear* In the selective integration 
technique, by using lower order integration for TSa y in 
Equation (8.3), the shear strain reduces to = *3 + > 
thereby retaining the pure bending behavior. Also for the 
Wilson's element in rectangular geometry, the 
contributions from and ^2 support the pure bending 
behavior exactly. Furthermore, since unique solution 
exists for the governing equations of elasticity, the 
stiffness matrix must be identical for all these cases. 
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8 . 2 8-NODE QUADRATIC PLANE ELEMENTS 



ISOPARAMETRIC ASSUMED DISPLACEMENT ELEMENT 


For the 8-node plane element. 


y. 


i -o 


(?I3) 


= I 


) z 3 4 


i <'-!*)(•♦ Ill) i * s, 7 

i 7*H' + J;f) i= *,S 


The displacement modes are 

Ui ^. + °<ar *3*) + ^ + «5 J* 

V = * *<, 3 ?* 


(y i4) 

1*1* “ul l 1 


The strain modes, 

e x -- «z ^ y* 

€y = «„ + <* lz * + Z* /+ y + X* +z^, 6 xy 


(f.is-) 
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*xy - + *'<>) + Cd+tZoc^) * (<*, x + z <=0 V 

* {Zoi s t- 2* {S )xy * * 7 x l + * u y* 

From here on the purpose for the isolation of the strain 
modes is to systematically construct the hybrid elements. 

HYBRID ELEMENT 

Following the procedure described previously for the 
linear elements f return to Equation (8.15) and eliminate 
artificial coupling and redundancy. 

0-* * /?, + fax +jS l0 xy 4 * /3 y x 

. * O'O 

<T y * fa + fa* + fay + fa *y + /*I4 x 

o-^ = ♦ fit x */S, y v* xy + N ** r yI 

The equilibrium equation yields following constraints on 
the ft 9 s 

A + A " ° 

* A = 0 (S./7) 

Al s 0 
#0 + ^.6 =° 

/£|| + = 0 

Applying the constraints and shifting the coefficient 
numbers yields following stress modes. 
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°* - /*■ * P+* + 

= /Sz **■ /ff 7 y + /£)*y 

oi y -- /S 3 -^/ -/s 7 x-i^ # y 

Since a minimum of 1 3 /S f s are required to obtain unique 
solution and thus eliminate the zero-energy deformation 


modes, two extra 

stress modes have 

to 

be included in 

the 

Equation 

(8.18). 

Furthermore , 

for 

the general purpose 

element , 

the two 

extra modes serve 

no 

purpose other 

than 

the elimination 

of zero-energy 

deformation modes* 

As 


before, the two modes are classified as the "nonsense" 
stress modes. With this in mind, all constraints imposed 


by the governing 

equations can be 

relaxed 

for the 

"nonsense" stress modes without loss of 

accuracy • 


Two possible 

candidates, in view of 

numer ical 

integration and 

application of 

uncoupled 

stress 

formulation, are 




c; -- 

/*.«. x y i 


(?. lSo) 


= *v 



These modes interact with E^x and £y*y 

to form 

non-zero 


energy* For more in depth coverage refer to the 
reference [10]. 


48 



8.3 8-NODE LINEAR SOLID ELEMENTS 



ISOPARAMETRIC ASSUMED DISPLACEME NT ELEMENT 
For the 8-node solid element, 

u. = a/ £ 

The displacement modes are 

u: + + «i'IT + ' t 7lT'- «t J7T 

v s + + *»s TV +0< ifJ1 +0< mT 7f 

w-s o < l7 ♦ *„ 7 + **0 5 4 >7 1 -°‘«3r + ^ f T 7 j 


The strain modes, 

£ x = ** + *s y + * 7 ^ *• y 2 

6y = +• °<ljX * *14 * + 

£ 2 = cx a<> + « ai 7 + <*^X + <**4 *y 


(S’. 20) 


49 



Y xy = (* 3 + *.o} + + «| 3 y + (*< + *, 5 )* + °^Y * Z *<*„ ys 

Ty^ = (c / (1 + <X IS ") + (ctT^ + X •+■ o( |4> y + o< zl l + <* u xy + 0(^x2 

r (*4 tf W) ‘*' 0< 7* + (*<+°(ii)y + ** 3 2 -f <*,*y + «* 4 y2 


HYBRID ELEMENT 

First eliminate the artificial coupling 
equation (8,20). 

<r x « fa + fay + fa a * + A* v* 

Cy ~ Pi + A* + fat z * A* Xz 

~ fa 4 A ^ ■*■ ^i« / + A*' x y 

*^y = fa * $3* + fit y + A* xz + + fa % ? 

°yi * A + A+ ? + fat* * /4. x y + A/ Xi ^ A*» x 

5ir^ * fas*+fa 2 +A4 x y - At y* + ^Jo/ 

From equilibrium, 

A< + A» -° 

/a + A =o 
As- * 0 

A* = Aa s A 4 = A* = A* "A 7 * 0 

Applying these constraints, 


in the 


(?.<2l) 


(*.**) 
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iS. Z$) 


o; = fy + /S u yi 

®V = /**. t / s »^^ 1 2 f / S , 7 

= ^3 + Ai x + *y 


°*v = 

Vl3 l 

4 /w 

4 


®y i z /* s 



1 

« 

4 ^a«y 


4 ^‘5 y 

" ^9 2 




From the three dimensionality of the element 

construction, an additional step must be performed before 
arriving at the final stress modes* By comparing the 

Equation (8.23) with the Equation (8.20), the modes 
^ , and j3 zi above are not represented as a possible 
deformation under the isoparametric formulation. Thus, 
these three additional fi's must be set to zero. Note that 
through the process of eliminating the artificial coupling 
in the Equation (8.21), additional possible modes under 
equilibrium are introduced. These ^ ' s do not contribute 
in accordance with the deformation modes and thus should 
be left out. 

The resulting stress modes are 

- P ' 4 Pi y fro* + fin'll 

- Px + Pi * + 2 . + ? n X 1 

^2 = Pi r pi x 4 p,z y + /*iy xy ( *•*+) 
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a + 

^13 ^ 

Hi 

II 

/** ♦ 

Ah * 

il 

lo* 

A + 

/*«y 


Recall from the Section 6 and the Figure 
and are "nonsense" stress modes 

equations can be relaxed for these terms 


6.1 that fl t , , 

and the governing 


UNCOUPLED STRESS HYBRID ELEMENT 

-""•a- 

In order to achieve fully equivalent element as the 
previous hybrid element with stress assumption in Equation 
(8.24) following uncoupled stress assumption must be used 
for computationally efficient element with isotropic 
material property* 

c; - * /Ay -+/?4-2 + *y + {*& yz * fa *2 

" P * + fa x + /^i* y ^ fa a + A* x y * A* + / rf /4x« 

= As + + /«7 y + 2 + A i *y ** A* y* x? 

cs; Y * A* + fai b 

°V* - A-4 *» As x 
<5 * t ■ Ai + A 7 y 
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Proper ^6' s can be eliminated or coupled directly by 
choosing the R matrix according to 

? T /! * o <>• 20 


The order of R matrix for the above stress assumption is 
number of ^8 1 s by 9. Thus additional inversion of (9*9) 
matrix ( R T R ' R ) will be required* 

mm ** 

In the natural coordinate system of the element, the 
stress assumption corresponding to Equation (8*25) is 


»■« * A ’■ M * hi +a t *fsm f A77 r A rr 

°V = A *• fa 7 * A 1 ) V" 7 + A* T? ’A 7 j * 

S = A * A 5 - A 1 + A S * A ?7 *A*7I + A< ?r 

1 A* ** A* J 
®y* " A« * As 7 
°5» = A< + A’ 7 


with u ^ as 

u x = A >• A 7 + ^ s T 
A = A ♦ A 7 ’ A f 
A = > 7 + >,7 *• A 7 


(3.27) 


Note that if the element is rectangular, above stress 
assumption reduces equivalently to Equation (8*33) since 
the Jacobian is constant. 

A large reduction in computional cost can be induced 
by recognizing that the quadratic terms are the "nonsense" 
stress modes and serve only to suppress the zero-energy 
deformation modes. With this in mind a resulting stress 
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assumption arrived is given below 


= h * fa hl'lS + TV 

e-y * P‘ * A t t fa 7 t *■ /?.. f ti+k+jv 

°i =A» * V'» 7 t / St »T + f J7*9f + jy) 


% -- fa ' f>n y 

°Vi s /i? + P'1 y 
~ f*o + P 1 ' 7 

With 



Thus the order of (rTi^R) is 


(?. zs ) 


(ff. z<\) 


reduced to (3x3) and the 


order of other matrices necessary for the stiffness matrix 
generation have been correspondingly reduced. 


TThw 

Recall that for the anisotropic material 
order for 

O' - C £ 


law, in 


(2.S) 


the constraints on cr 


and £ are 



( 7 . 29) 


P 


i 


ft 


p 


3 




(7. 
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A promising choice under this constraint for the 


stress 


a ssumpt ion i s 


°; = 

A 

+ M 

- A 7 

-At 

°y = 

A 

"A! 

-A 7 

-A T 

* 

A 

4 A I 

’A* 7 

- at 

V 

A* 

+ A»T 

-A»7 

-AT 

C Yt = 

A, 

+ AO 

-A>1 

-A-T 


: 

+ A*7 

- A" 7 

-AiT 


+ 1**0 (si + TO 

+fi*(Ji+iT + 7r) 

(?.30) 

*faoL SI + 11 + n) 

+ p 2 s(Sl+'JJ+Jf) 

+ fizo ^7*7r + jO 


However, since the stress distribution in Equation (8.30) 
may be too rigid, additional relaxation of the governing 
equations may be necessary. 
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8 . 4 20-NQDE QUADRATIC SOLID ELEMENTS 



ISOPARAMETRIC ASSUMED DISPLACEMENT ELEMENT 


For 20-node solid elements. 


U * 


(?. 31) 




4 ( 1 - 7 1 ) C^J L J) 

i ('-x'Hi+JiV O7.7) 


1 * 10 , 11 , 1 4,»* 

i s l Z iV%*“ 


The displacement modes are 
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The breakdown of the corresponding strain modes are 
included in the Appendix A. 
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HYBRID ELEMENT 


The procedure for the 20-node solid hybrid element 
development is similar to 8-node solid hybrid element. 


The stress modes 

are given 

in the 

Appendix A. 

In 

the 

final resulting 

assumed 

stress 

mode s , 

Table 

A 4 , 

the 

"nonsense" modes 

include all cubic 

stress 

terms • 

This 

i s 


consistent with the linear 8-node solid element where all 
quadratic terms are the "nonsense" modes. Furthermore, 
this is mathematically consistent since with quadratic 
displacement assumption, cubic stress modes does not 
contribute to general convergence of the element. Recall 
that both the stress and the displacement approximations 
must be improved properly and simultaneously. Since the 
accuracy of stresses from the displacements converge, at 
best, in a quadratic order (Optimal Gauss Points [ 1 7] ) , the 
cubic terms are only used for suppressing the zero-energy 
deformation modes. 


57 



9. NUMERICAL EXAMPLES 


The element nomenclature is given in the Appendix B 
for quick reference. In order to numerically support the 
discussion in the previous sections, various samples of 
stress assumptions have been implemented. However, bear 
in mind that most of the numerical examples given are 
carried out to provide a foundation for the inductive 
process used in the previous sections. 

All calculations for the numerical examples are done 
under double precision using Digital Vax 11/780 computer. 
The Gaussian quadrature order used for the numerical 
integration correspond to exact integration for each 
element type . 
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9 . 1 CANTILEVER BEAM USING SINGLE 4-NODE PLANE 


STRESS ELEMENT 

Refer to the Figure 9.1 for the problem description 
and geometry. The cantilever beam problem has been chosen 
to depict the bending behavior. Since the purpose for 
this example is to illustrate the nature of the 

interaction of the stress modes, only single element mesh 
is needed. The Poisson's ratio has been set to zero to 
isolate the modes. The analytical solution is based upon 
the Bernoulli beam theory. The tip deflection and the 
maximum normal stress, 0 ^ , for the five different stress 
assumption is given in the Table 9.1. 

TABLE 9.1. Cantilever beam using Single 4-node element. 


ELEMENT 



RP4A 

5 A's 

.015 

3000 

RP4B 

5A's 

.0012 

0 

RP4C 

7P's 

.0011 

226.3 

RP4D 

9 P's 

.0011 

222.2 

RP4E 

SP's 

.0012 

0 

DP4 (7T P ) 

.0011 

222.2 

Analytical 

.015 

3000 


The non-zero /5's for each stress assumptions for the 
cantilever beam bending are listed below. 
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RP4A 


=/*, v*y 


RP4B 


RP4C 


RP4D 


RP4E 





p*+r$ y 
l <W» -fis* 


Note that only in the Element RP4A the pure bending 
behavior is exactly modeled. Thus even in a complex mesh 
model, the transmitted bending load to each element is 
successfully modeled by the Element RP4A. In any complex 
loading problem the bending load will be present for many 
e 1 ement s • 

The reason for the choice of other stress assumptions 
are as follows: 

RP4B - 5 p case with zero-energy deformatfion mode 
suppr e s sed • 

RP4C - Complete linear stress assumption with 
equilibrium condition imposed. 
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RP4D - Complete linear stress assumption. 

RP4E - Alternate case with equilibrium satisfied 

As indicated by the results, if the corresponding stress 
modes are available, the interaction between the stress 
and strain modes will reintroduce the artificial coupling. 
By using the minimum required number of ^*s, 5 in this 
case, and choosing the stress modes in full recognition of 
the strain modes, above interaction can be avoided. 

The nature of the stiffness matrices for a square 
element can be summarized by the examination of the trace. 

TABLE 9.2. Trace of the stiffness matrix. 


ELEMENT 

Trace = 2! Eigenvalues 

RP4A 

.3634 x10* 

RP4B 

.3223 

RP4C 

.3727 

RP4D 

.3956 

RP4E 

.3152 

DP 4 (-n>) 

.4615 


The Element RP4E is most flexible element in sum, yet this 
flexibility does not extend to the bending behavior. 
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9.2 CANTILEVER BEAM USING SINGLE 8-NODE PLANE 


STRESS ELEMENT 

Identical problem shown in the Figure 9.1 is solved 
using single 8-node plane stress element. Recall that 

during the process of the derivation of the hybrid 8-node 

plane element/ the strain modes from the displacement 

modes/ Equation (8.15) , has been shown to contain the 

bending nodes, °^4 and Thus any stress assumption with 

bending modes must give exact solution. The tip 

deflection results are given in the Table 9.3. 

TABLE 9.3. Tip deflection - Cantilever Beam Problem. 


ELEMENT 

^TIP 

DP 8 

.0137 

RP8A 

.015 

RP8B 

.015 

RP8C 

.015 

RP8D 

.015 

RP8E 

.015 

RP8F 

.015 

Analytical 

.015 


The assumed displacement element, DP8, did not give the 
exact solution due to the interaction of ^ and ^ tern in 
the Equation (8.15). As a point of interest, the trace of 
the stiffness matrix from the Elements RP8E and RP8F is 


identical 
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9.3 CURVED CANTILEVER BEAM USING 8-NODE PLANE 


STRESS ELEMENT 


The curved cantilever beam problem f Figure 9.2, has 
been motivated by Spil ker , Ma s ker i # and Kania[16] • In the 
article, to achieve invariance under coordinate rotation 
the stress modes are expanded to full cubic. In order to 
further reduce the number of ’s, the compatibility 
condition is imposed on the stress. In Appendix B, this 
element is named RP8D • The Elements DP8 , RP8A, RP8B , and 
RP8C are used for comparison. 

The tip deflection results are shown in the Figure 
9.3. The analytical solution has been obtained from 
Timoshenko and Goodier[18] . Note that this solution is 
approximate and thus the percent error only has an 
approximate meaning • 

The stress results for the five element mesh are 
plotted on the Figure 9.4. The results are obtained at 
various angles for r* 1 1 • 58 . Specifically# these points 
correspond to the optimal stress points. Since the stress 
distributions obtained from the Elements RP8A , RP8B , and 
RP8C are close, only the RP8A stress distribution is 


shown • 

Overall results indicate 
converges much faster then 
formulation in the curved beam 
the Hybrid elements there is no 


that the Hybrid formulation 
the assumed displacement 
analysis. However, between 
consistant way to evaluate 
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which element has the best result* To eliminate the 
occillation of the Element RP8 A stress distribution two 
additional analyis should be made. 

1) Use local coordinate for RP8A, RP8B, RP8C Elements. 

2) Develop an element using the natural coordinate 
system for the stress modes. 
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9 • 4 CANTILEVER BEAM USING 8-NODE SOLID ELEMENT 


The cantilevered beam problem, Figure 9.5, is solved 
using several mesh arrangements. Figure 9.6. A moment. 
Load Case I, and a shear. Load Case II# loadings are 
considered. The analytical solution is based upon the 
Bernoulli beam theory. The tip deflection results are 
given in the Tables 9.4 and 9.5. 

In a rectangular solid mesh configuration, the tip 
deflection results demonstrate that the Elements RUS8A and 
RUS8B contains the pure bending modes induced by the 
moment couple in the Load Case I. As expected, the 
results from the Elements RUS8D and RUS8E are similar to 
the result from the assumed displacement element DS8. 
This is a further evidence that inclusion of unnecessary 
P's simply reduces the hybrid element behavior similar to 
the assumed displacement element. 

In order to assess the rate of degredation of 
accuracy as the elements are skewed, the cantilevered beam 
problem is repeated by steadily increasing the distortion* 
The result from the distortion sensitivity analysis of the 
hybrid element, RUS8A, is plotted on the Figure 9.7. As 
shown, the Element RUS8A has a high degredation rate 
initially and then effectively reduces to zero after b/a = 
2. An interesting point to consider is that the value of 
the tip deflection beyond b/a = 2 is approximately the 
same as the result obtained from the assumed displacement 
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TABLE 9.4. Tip Deflection 


LOAD CASE I 


ELEMENT 

MESH1 

MESH2 

MESH3 

MESH4 

MESH5 

DS8 

9.00 

27.78 

20.27 

23.52 

18.81 

RUS8A 

100 

100 

51.1 

46.0 

25.8 

RUS8B 

100 

100 

51.1 

43.5 

25.8 

RUS8D 

9.26 

30.2 

22.5 

23.0 

16.6 

RUS8E 

9.26 

30.2 

22.5 

23.3 

16.7 

Analytical 

100 

100 

100 

100 

100 


LOAD CASE II 


ELEMENT 

MESH1 

MESH2 

MESH3 

MESH4 

MESH5 

DS8 

9.26 

28.50 

22.81 

24.08 

20.77 

RUS8A 

77.5 

96.0 

58.9 

55.2 

40.1 

RUS8B 

77.5 

96.0 

58.9 

53.3 

40.1 

RUS8D 

9.44 

30.8 

23.8 

25.4 

20.8 

RUS8E 

9.44 

30.8 

23.8 

25.7 

20.9 

Analytical 

100 

100 

100 

100 

100 
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element, DS 8 . Basically, the distortion reintroduces the 
artificial coupling into the stress modes* This 
observation provides a clue to develop a more effective 
element with less sensitivity to distortion* The tip 
deflection result for b/a * 99 is 34*9 which further 
supports that the accuracy of the hybrid element will be 
equal or greater then the assumed displacement element no 
matter how large the distortion becomes for the 
cantilevered beam problem. 

Backtracking for a moment, re-examine the Table 9.4 
and 9.5 comparing the Elements RUS8A and RUS8B. Recall 
that these two elements employ similar assumed stress 
modes. Only difference being that one uses the xyz 
coordinates and the other the natural coordinates. The 
comparison of these elements indicate further modification 
is necessary to achieve distortion insensitive element 
then just expressing the stress modes in the natural 
coordinate system. Above remark seems reasonable since 
the mapping a linear mode from xyz coordinate to the 
natural coordinate system yields also a linear mode. As a 
point for future research, the mechanics of distortion 
should be studied in view of the coupling between the 
stress modes • 

In any structural problems, a reference coordinate 
system must be defined to locate each point in the 
volumetric space. Thus next step of analysis studies the 
solution response of the cantilever beam problem under 
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concentrated tip load as the reference coordinates are 
rotated. Since the stress assumption for the Elements 
RUS8A and R0S8B are not a complete polynomial r the 
invariance condition is not automatically satisfied. 

Before presenting the results, note that an alternate way 
to establish the reference coordinate invariance is to use 
a local coordinate system for each element. 

The results from the rotation of the reference frame 

is shown in the Figure 9.8. The Element RUS8A using the 

Cartesian coordinate system for the stress inodes contains 

A o 

zero-energy deformation modes at w =* 45 • This is 

verified by eigenvalue analysis of the stiffness matrix 
generated at this angle. Furthermore, as the angle 

approaches the value of 45° the element becomes more and 
more flexible till it becomes unstable at © =* 45 . 
However, by the use of the natural coordinate system, 
Element RUS8B , eliminates this unstable mode. Notice that 
the variation of the result remains negligible for 0< 
30° . 

Before the reader becomes too astounded by above 

results, several comforting observations are in order. In 
most engineering problems, the reference coordinate chosen 
coincide with the physical structural geometry which 
eliminates the possibility of the complete structural 

instability. To clarify, due to the boundary conditions 
and the assembly with stable elements, the total structure 
will be stable even when part of the mesh is parallel with 
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possibility 


45° line. Another strategy to eliminate the 
of the unstable mode is partially distorting the element 
in the same plane of the expected angle 0 . Following 
result demonstrates this numerically. 


. o 

TABLE 9.5. Tip deflection, Load Case II at <p *45 . 


ELEMENT 

MESH2 

MESH3 

MESH4 

RUS8A 

-.599X10 13 

-.240*10 5 

132 

RUS8B 

159 

143 

129 


The angle lies in the x - y plane, Figure 9*5. Since the 
Mesh 4 is distorted in the x-y plane, the unstable mode is 
partially dampened. Above analysis also holds for the 
reduced integration technique. Same approach can be used 
to employ elements with zero-energy deformation modes. 

As a final note, for the quadratic elements, since 
the stress assumption is complete to the linear order, no 
instability will arise for the bending problem. Overall, 
the remedy to introduce invariance with respect to the 
reference frame should be made by using the local 
coordinate system. The problem of arbitrary reference 
coordinate system should not enter into the element 
development. Combat coordinates with coordinates. 
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9*5 CIRCULAR HOLE IN AN INFINITE STRIP USING 8-NODE 


SOLID ELEMENT 

The model for the circular hole in an infinite strip 
is shown in the Figure 9.9. The Figure 9.10 provide a 
plane view of the four configurations of mesh using 8-node 
solid elements. Although, due to the curved geometry of 
the hole, higher order element should be used, this 
problem nicely demonstrates the limitations of the hybrid 
elements developed for general purpose applications. 
Furthermore, this problem establishes the groundwork for 
the discussion of the special purpose elements presented 
in the next section. 

The result from the displacement convergence study is 
given in Figure 9.11. The stress distribution along x * 0 
for the four mesh configurations are provided in Figures 
9.12 to 9.15. The elements used in the analysis are the 
Elements DS8 and RUS8A. Recall that RUS8A can model pure 
bending exactly as demonstrated in the subsection 9.4. 
Since the results from the Element RUS8B with stress 
assumption in the natural coordinate system are very close 
to the results from the Element RUS8A, these results are 
omitted • 

In the overall sense the difference between the 


results 

from 

Element DS8 and 

RUS8A are not 

overly 

dramatic . 

In 

the development 

o f 

the 8-node 

sol id 

elements , 

the 

assumed stress modes 

clearly show 

that the 
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over 


the 


assumed 


advantage of the hybrid 8-node 

displacement 8-node is the bending behavior. In any other 
type of mode excitation, the convergence of the 8-node 
solids using either element will be similar. However, 
note that this characteristic is purposely imposed in 
order to construct a general purpose element. 

In the stress distribution obtained, the results from 
the assumed displacement element are actually little more 
accurate then the result from the hybrid element. The 
only way this can be explained is that the artificially 
coupled terms inhence the accuracy for the class of 
problems exemplified by the circular hole problem. From a 
logical extension to the above conclusion, introducing 
more descriptive modes will increase the accuracy for the 
corresponding class of problems, i.e. special purpose. 
Further discussion is reserved for the next section. 
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9.6 HOLLOW SPHERE UNDER TEMPERATURE DISTRIBUTION 
USING 20-NQDE SOLID ELEMENT 

The problem chosen to evaluate the 20-node solid 
elements is the hollow sphere under temperature 
distribution, Figure 9.16. Six element mesh is adequate 
to study the relaxation and symmetry condition on the 
cubic terms. Six diffenent stress assumptions are used 
for direct comparison. Element used in this analysis are 
listed below with comment for their choice. 

ASSUMED DISPLACEMENT ELEMENT 

DS 2 0 

HYBRID ELEMENTS 

RS20A - Equilibrium relaxed for 

all cubic stress modes. 

Symme try Maintained. 

? % 

RS20B - Equilibrium relaxed for only <7 X = x , (Ty *y t 
and CT t » z* term. 

Symmetry Maintained. 

RS20C - Equilibrium imposed in unsymmetric fashion. 

RS20D - Equilibrium relaxed for both quadratic and 
cubic stress modes. 

Symmetry Maintained. 
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RS20E - Alternate stress modes used to complete the 
stress assumption, 57 *s. 

Equilibrium relaxed for 

all cubic stress modes. 

Symmetry Maintained. 

RS20F - Also an alternate form, 54 1 s. 

Equilibrium relaxed for 

all cubic stress modes. 

Symmetry Maintained. 

The analytical solution for the hollow sphere problem 
is obtained from Timoshenko and Goodier[18]. A comparison 
of the radial displacement distribution obtained by 
Elements DS20 and RS20A is shown in Figure 9.17. The 
results are both very accurate when compared with the 
analytical solution. Again for the tangental stress 
distribution result, Figure 9.18, both elements provide 
accurate solutions. This is expected since the rate of 
curvature change is slow for the tangental stress 
distribution. Thus the excitation of the higher order 
stress modes are correspondingly small. 

The radial stress distribution provide high enough 
rate of curvature change to excite the higher order stress 
modes in order to distinguish each element performance. 
The Figure 9.19 compares the radial stress distributions 
obtained from the assumed displacement element, DS20, and 
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the hybrid element, RS2 0 A . In general, a significant 
disparity of results, exemplified in Figure 9.19, between 
/ \p and TTh element will be observed for any problem 
with stress distribution that has a high rate of curvature 
change. In another words, when the quadratic stress modes 
are excited disparity between the two formulation will 
arise. Above remark is a simple extension of bending 
problem, linear mode excitation, using 8-node solid 
element. 

The radial stress distribution comparison between 
hybrid elements with various other stress assumptions are 
shown in Figures 9.20 to 9.23. The results indicate that 
the equilibrium can be relaxed for all cubic stress modes 
and the symmetry condition should be maintained. Also, 
the elements RS20A, RS20E, and RS20F gave almost identical 
stress distributions. Overall, the hollow sphere problem 
numerically supports the discussions on the previous 
sections • 
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9.7 CANTILEVERED BEAM USING 20-NQDE SOLID _ ELEMENTjj 


In a recent article by Spilker and Singh[19] , a 
hybrid element with complete cubic assumed stress 
distribution satisfying both equilibrium and compatibility 
conditions is presented. This element is named RS20G 
(refer to Appendix B). An example given is a cantilevered 
beam problem with distributed end shear loading. The 
normal and shear stress distribution results are given in 
Figures 9.24 and 9.25, respectively. The results from all 
three elements, DS20, RS 2 0 A , and RS20G, are comparable as 
expected from previous beam bending analysis. 

The purpose for constructing the Element RS20G, by 
Spilker and Singh, is to implement element invariance with 
respect to the reference coordinate. The invariance is 
accomplished by expanding the stress into full cubic 
distribution. The stress compatibility condition had to 
be used to reduce the number of P • s to 69. This is 
compared to 54 P'S used for the Element RS20A. 

Under a close examination, several severe limitations 
on the Element RS20G restrict its applicability. First, 
the element is limited to isotropic material only. Also, 
excessive computational cost prohibit practical 
application. Uncoupled stress formulation cannot be used 
to reduce cost since 51 constraints are necessary to 
reduce full cubic, 120 P's, to 69 P's. A lesson learned 
from this example is to use local coordinate system to 
achieve element invariance. 
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9 . 8 HU-WASH IZ U ELEMENT EVALUATION 


In order to numerically compare the difference in the 
solution between t and when the stress-strain 
relation is relaxed, two previous problems, given in the 
subsections 9*6 and 9.7, are solved using similar P 
matrix. The elements used are 20-node solid elements 
RS20A and HS20A (Appendix B). 

For the cantilevered beam problem, the tip deflection 
results are given in Table 9.6. 


TABLE 9.6. Tip deflection of cantilevered beam. 


# of Elements 
in the mesh 

1T ? , DS20 

, RS20A 

HS20A 

1 

-. 956 * 10 '* 

-1.05*10** 

-1.05*10* 

2 

-1.08 

-1.09 

-1.09 

4 

-1.10 

-1.11 

-1.11 


S analytical = ~ 1 - 132 * 10 ~* (Transverse shear included) 


Also the shear stress results 
identical up to 3 digits, 
graphically. The reason the 
HS20A are similar is that 
has a linear distribution in 
constant distribution in the 
of the stress assumptions for 


between RS20A and HS2 0 A are 
Figure 9.26 demonstrates this 
results between RS20A and 
the cantilevered beam problem 


the 

normal 

stress and a 

shear 

stress • 

In comparison 

the se 

two elements. Appendix 
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B, showes that for the constant and linear terms, the 
equivalence conditions are identically satisfied. 

The previous discussion on the Hu - Wa shi zu formulation 
asserted that the str e ss- stra in relation can only be 
relaxed for the cubic "nonsense" stress modes. To 
illustrate the consequence when the stress*strain relation 
is also relaxed for the quadratic modes, the hollow sphere 
problem is solved using Element HS20A. Recall that the 
condition required to satisfy the stress-strain relation 
constrains the 1 s to be equal. The Figure 9,27 clearly 
demonstrates the difference between the two formulations. 
Since the stress-strain relation for the quadratic terms 
for Element HS20A are not satisfied, the result is poor. 
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9 . 9 CPU TIME FOR STIFFNESS MATRIX GENERATION OF 


8-NODE SOLID ELEMENTS 


To justify the use of the hybrid elements for 
practical industrial application, normalized CPU time for 
stiffness matrix generation of 8-node solid elements are 
provided in Table 9.7. 


TABLE 9.7. Normalized CPU time for stiffness matrix 
generation of 8-node solid elements. 


DS8 - Assumed Displacement Method 1 

RS8A - Original Hybrid Stress Method 1.62 

RUS8A - Uncoupled Stress Method 0.96 

RUS8C - Uncoupled Stress Method 0.60 


Note that all three hybrid elements, RS8A, RUS8A, and 
RUS8C, provide exact pure bending behavior. Thus by using 
the flexibility allowed in the uncoupled stress method, 
the economic roadblock on the hybrid elements can be 
easily bypassed. 
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10. SPECIAL PURPOSE ELEMENTS 


The special purpose elements, as differentiated from 
the general purpose elements, are tailored for a 
restricted class of problems. By restriction, the 
convergence using the special purpose elements can be 
significantly increased. The aim of the present 
discussion is to direct future research in the special 
purpose elements. 

In the previous sections, much emphasis is placed on 
the examination of both displacement and stress modes in 
conjunction. The resulting combinations that satisfy the 
governing equations establish the accuracy of each 
element. If the exact mode is included in the assumed 
modes, only a single element is necessary to obtain the 
exact solution. If otherwise, the convergence depends on 
the ability of a linear combination of modes to 
approximate the solution. 

Inevitably, the next generation of elements will 
involve a joint effort of theoretical solid mechanics and 
finite element strategy. Expanding, as Wilson attempted 
to purposely insert bending modes, the modes obtained from 
analytical means will be implemented into the assumed 
modes. Observe that the displacement and stress modes 
obtained analytically for a sepecific problem satisfy the 
governing equations a priori. Furthermore, for a similar 
class of problems, this mode will be highly excited. 
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Overall, above technique simulates the Eigenmode expansion 
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nsion used for many years to solve 
ms for the finite element method, 
lity is established and thus the next 
feasibility. To begin with, the use of 
ation provides a convenient method to 
sired modes. By examining, step by step, 
of hybrid element, the only roadblock 
numerical integration of the elementary 
as the trigonometric and exponential 
rute force solution is integrating each 
y. A tremendous amount of necessary 
be done using an algebraic manipulator, 
ill be found as the research evolve in 
f special purpose elements. 
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1 1 . CONCLUSION 


The key to unlocking the mystery of the Hybrid stress 
formulation is to understand the interaction of the 
assumed displacement and stress distribution. Using the 
isoparametric method as the basis for the assumed 
displacement, a consistent stress distribution can be 
obtained. In practical engineering applications, each 
element is employed to cover a "finite" domain. 
Therefore, the mathematical proofs on convergence, derived 
on the assumption of the limiting process, have only 
restricted practicality. The bridge to fill the gap 
between mathematical abstraction and practical reality of 
finite element method is the application of the physics 
involved in continuum mechanics. In this context, all of 
the painstaking development of analytical solution can be 
applied in the finite element method. Recognize that the 
finite element method is a numerical method operating on 
the parameters provided through the assumed variable 
distributions. 

In the discussion of general purpose elements, the 
Statement of Equivalence convey that if the assumed 
displacement, obtained through isoparametric formulation, 
is complete to order P, then the Hybrid and the assumed 
displacement element stress accuracy is equivalent to 
order (P-1). The actual superiority of the hybrid element 
will be visible for the problems requiring the stress 
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modes represented in the order P. This is demonstrated by 
constant stress problem verses cantilever beam problem for 
8-node solid element (P*1). 

In sum, the most attractive attribute of the hybrid 
stress formulation has not been yet fully exploited. The 
flexibility to easily include assumed modes obtained 
through analytical methods differentiates hybrid elements 
from the assumed displacement elements. Always remember 
that flexibility is an attribute when applied with 
understanding . 
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FIGURE 6.1 Possible deflection mode for 8-node solid 
requiring a "nonsense" stress term (in this 
case0i= xy) to prevent such a zero-energy 
deformation mode [15]. 
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P 11=0 

u=v=0 

P = 1000 
E = 10 7 
0= o 

FIGURE 9.1 Cantilever beam using single 4-node plane 
stress element. 
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FIGURE 



E = 10 7 r a = 12 

4 = .3 r. = 10 

t 


F 1 = 270.3 

~ 56.3 Total Moment = 600 
F 3 = 326.6 


MESH 

1 

= 16 DOF 

( 1 element) 

MESH 

2 

= 26 

( 2 elements) 

MESH 

3 

= 36 

( 3 elements) 

MESH 

4 

= 46 

( 4 elements) 

MESH 

5 

= 56 

( 5 elements) 


9.2 Curved cantilevered beam using 8-node plane 
stress elements. 


85 




80333 4 'N0I1031J3a dll 


o 

Q 

LU 

LU 

C£ 

Lu 



86 


FIGURE 9.3 Tip deflection for curved beam. 
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FIGURE 9.4 Tangental stress at optimal stress points, r=11.34, 
for 5 element mesh curved beam. 
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FIGURE 9.5 Cantilever beam using 8-node solid elements. 
(Single element mesh - MESH 1) 
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FIGURE 9.8 Deflection of cantilever beam under concentrated tip 
load as the reference coordinates are rotated. 
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FIGURE 9.9 Circular hole in an infinite strip. 
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FIGURE 9.13 Stress distribution along x=0, MESH 2. 
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FIGURE 9.14 Stress distribution along x=0, MESH 3. 
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FIGURE 9.15 Stress distribution along x=0, MESH 4. 
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FIGURE 9.16 Hollow sphere under temperature distribution 
using 20-node solid elements. 
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FIGURE 9.17 Radial displacement distribution along the x-axis. 



10E3 



101 


FIGURE 9.18 
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FIGURE 9.19 Radial stress distribution. 
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FIGURE 9.20 Radial stress distribution. 
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FIGURE 9.21 Radial stress distribution. 
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FIGURE 9.22 Radial stress distribution. 
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FIGURE 9.23 Radial stress distribution. 
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FIGURE 9.24 Cantilevered beam using 20-node solid elements 
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FIGURE 9.26 Cantilevered beam using 20-node solid elements. 
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FIGURE 9.27 Hollow sphere under temperature distribution. 



APPENDIX A 


TABLE AT. Strain modes from assumed displacement modes. 
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TABLE A2. After elimination of the Artificial Coupling. 
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TABLE A3. Equilibrium applied. 



114 



















































































When the Table A4 is compared to the Table Al, the 
stress modes fi So * ftsz* and fss are not present. Thus these 
modes are due to the process of decoupling the stress modes. 
Furthermore , the presence of and ^ are redundant 

and should be removed. In the final form of the stress 
assumption i Table A4, the equilibrium condition is relaxed 
for the cubic terms. Also for convenience, the numbering 
scheme has been changed. The x^, y 3 , and z 3 terms are 
chosen to eliminate the additional zero-energy deformation 
mode s . 
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TABLE A4. Final form with equilibrium for cubic terms relaxed. 
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APPENDIX B 


ELEMENT NOMENCLATURE 


ASSUMED DISPLACEMENT FORMULATION 


D - Assumed Displacement 
P - Plane Stress 
S - Solid 

# - Number of Nodes 

DP4 - 4-Node Plane Stress Element 
DP8 - 8-Node Plane Stress Element 
DS8 - 8-Node Solid Element 
DS20 - 20-Node Solid Element 


HYBRID FORMULATION 

R - 

» - TThw 

U - Uncouple (Blank for Coupled) 

P - Plane Stress 
S - Solid 

# - Number of Nodes 

A,B,C,*.. Assumed Stress Version 
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4-Node Plane Stress Elements 


RP4A 


RP4B 


RP4C 


RP4D 


RP4E 


•/, * fay 

°V = fa + fa * 
* fa 

<*k * fa 

°V r A 





v*y V'X 

oy = ys s t fax ■* fa y 7 P 5 

<V p s - fax- fay 


v< * fa +fa* + fay 

fly = fa'* fax. +fa y 
o- xy = + fa* + fay 


°< = Z 3 1 + ftx 

Gy* fa + fay s fas 

% = fa ~fai -fax 


fa i*) 
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8-Node Plane Stress 
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Const. 
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* 

X*Y 




X 3 
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RP8C 


Elements 

Delete circled terms. 
(Equation 8.22) 

All terms included. 


Alternate version 
of RP8A. 
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RP8D - Complete Cubic 
Equilibrium and 
Compatibility impose' 

(Ref. [16]) 


RP8E - Stress assumption 
to only suppress kinematic 
modes without regards to 
equil ibrium. 
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RP8F - Complete quadratic, ignores equilibrium. 
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8-Node Solid Hybrid Elements 


* Equilibrium condition is fully satisfied unless 
otherwise indicated* 

RS8A 

o; * A T 

«V Mi * 

"k s A » 

Ay = A- 7 
®ya = A " r 
Aa 1 fi + 


RUS8A - 9 Constraints on p's from Equilibrium Condition. 

<r x = >5, +^x ^ 3 y + /? 5 xy + +/ 7 x?. 

°V = Z 3 ? 4 fa * 4 y 4 2 4 /^ia*y 4 ^li yz 4 /^4X2 

s Af 4 4 A»7 / 4 ^ 4 Xy 4 ^ y E 4 X2 

<£y * Aa 4 * 

4 /^as* 

/*a* 4 y 


* Note that RS8A and RUS8A are equivalent element with different 

programming algorythm introduced in the uncoupled stress formulation. 


Ar * A.* + A*y* 

At X + An * + A»7 X * 

/,* * i*ij +^.r x y 

At* 

A** 

Au y 
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RUS8B - 9 Constraints. 


A * PzV m ^ Pa V + h yy Pi 71 * /Mr 

s + /M + AO + / S " T * P**.V1 + P^l^faJY 

^'5 + Pu^ + Pn*)* Pw\ + fa y*7 +/^ao7r *£ i( TT 

®xy r * Pz* y 

fly* 5 £*4 + /S 2S y 

°« = + /i-7 ^ 


*RUS8B is equivalent to RUS8A for rectangular elements. 


RUS8C - 3 Constraints. 
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*RUS8C - Equilibrium relaxed for quadratic terms. 
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RUS8D - 9 Constraints. 


G x s fi, + A* +/S 3 y + 

°V = A + A * + A / 

^ A + A X + ^.y 
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RUS8E - 9 Constraints. 
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RS20G (Ref. [19]) 


<r„ = 0, + 0 2 X +0J y +P*Z + 1/2 {y 2 -x 2 )0 2 2 + l/2A[-y 2 “Ax 2 + (1 + A)z J ]0u 
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